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Abstract 

Algebraic topology (homology) is used to analyze the weakly turbulent state of spiral defect 
chaos in both laboratory experiments and numerical simulations of Rayleigh-Benard convection. 
The analysis reveals topological asymmetries that arise when non-Boussinesq effects are present. 
The asymmetries are found in different flow fields in the simulations and are robust to substantial 
alterations to flow visualization conditions in the experiment. However, the asymmetries are not 
observable using conventional statistical measures. These results suggest homology may provide 
a new and general approach for connecting spatio-temporal observations of chaotic or turbulent 
patterns to theoretical models. 
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I. INTRODUCTION 



Recent technical advances in experimental fluid mechanics now make it possible to mea- 
sure complex dynamical behavior with high resolution in space and time (|l|). Similarly, 
modern computational fluid dynamics methods permit modelling of complex chaotic and 
turbulent flows (|2j). The data sets produced by such experiments and simulations can be 
enormous; as a result, interpreting the results becomes a significant challenge. In particular, 
characterizing the geometric properties of complex spatio-temporal patterns in large data 
sets has been difficult because, to date, no general methodology has existed for extracting 
geometric signatures. 

Algebraic topology provides a tool for describing global geometric properties of structures. 
Devised by Poincare 3| for use in global nonlinear analysis, algebraic topology originally 
used as input analytically defined objects (e.g., level sets of differentiable functions) to pro- 
duce an output in the form of algebraic quantities which convey topological information 
about the input. In modern times, input objects can be expressed either as simplicial or 
cubical complexes. In fluid mechanics and in most fields of science and engineering, cubical 
representations often arise naturally in both experiments (raw image data represented as 
square pixels or cubic voxels) and simulations (fields computed at gridpoints on square or 
cubic lattices). A package of computer programs has been developed to perform compu- 
tations of algebraic topology (computational homology) on cubical complexes in arbitrary 



dimensions. This suite of tools, cal 
available for download via the Web ( 



ed CHomP (Computational Homology Project) is freely 

13). 

We report here the first use of computational homology to characterize data obtained 
from a laboratory experiment. We analyze Rayleigh-Benard convection in the state known 
as spiral defect chaos ([6]), which is widely considered a paradigm for the little understood 
phenomenon known as spatio-temporal chaos ([3]). In planform, patterns of spiral defect 
chaos, which are observed just above convective onset in low Prandtl number (~ 1) fluids, 
are composed of convection rolls deformed into numerous rotating spirals and riddled with 
dislocations, disclinations and grain boundaries. Spiral defect chaos has been quantitatively 
described by a wide variety of approaches, including structure factors, correlation lengths 
and times as well as wavenumber, spectral and spiral number distributions ({§(]; see also 
and references therein). 
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Thermal convection is frequently modelled using the Boussinesq approximation, which 
assumes the temperature dependence of the fluid properties can be neglected, except for the 
temperature-induced density difference in the buoyant force that drives the flow. However, 
non- Boussinesq effects can arise in flows both in the laboratory and in natural settings. At 
convective onset, the subcritical bifurcation to hexagonal patterns is a clear signature of 
non- Boussinesq effects. (Straight convection rolls arise at onset from a supercritical bifur- 
cation in Boussinesq Rayleigh-Benard convection.) Non-Boussinesq effects can be described 
quantitatively using perturbation theory near onset; in this regime they are characterized 
by parameter Q introduced by loj]. Values of Q >~ 1 indicate significant non-Boussinesq 
effects; Q = for Boussinesq convection. As changes in control parameter move the con- 
vective flow well away from onset, non-Boussinesq effects typically become more important 
and are more difficult to characterize theoretically and experimentally. 

II. EXPERIMENTAL PROCEDURE 

We measure convective flow in a horizontal layer of compressed (3.2 MPa absolute pres- 
sure) CO2 gas of depth d = 0.0690 ±0.0005 cm. The layer is bounded above by a 5-cm-thick 
sapphire window and below by 1-cm-thick gold plated aluminum mirror. The lateral walls 
are circular, formed out of an annular stack of filter paper sheets 3.80 cm in diameter. An 
electrical resistive heater is used to heat the bottom mirror to a temperature, T^, and the 
top window is cooled to a fixed temperature of T c = 21.20 ±0.02 °C by circulating chilled 
water. When the temperature difference, AT = T h — T c , across the gas exceeds a critical 
temperature difference, AT C = 3.8 ±0.1 °C, the onset of fluid motion occurs. The Prandtl 
number Pr is 0.97. In the experiments, the system control parameter, the reduced Rayleigh 
number e = (AT — AT C )/AT C , is increased above onset through a range where spiral defect 
chaos occurs. The characteristic time scale, the vertical diffusion time t v , is approximately 
2 seconds. 



The flows are visualized using the shadowgraph technique ([llj). Time series of shad- 
owgraph images [Figured] (a,b)] with a spatial resolution of 515 x 650 pixels are captured 
under computer control at a rate of 11 Hz using a 12-bit digital camera interfaced to a 
framegrabber. The images are prepared for analysis by subtracting a background image of 
the fluid below onset and using digital Fourier filtering to remove high wavenumber compo- 
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FIG. 1: Images of spiral defect chaos convection are shown from laboratory experiments (a) & 
(b) and numerical simulations (c) & (d). Shadowgraph images from the experiments illustrate 
the convective flows at (a) e = 1.0 and (b) e = 2.0. The mid-plane temperature field is shown at 
e = 1.4 for simulations carried out under (c) Boussinesq (Q = 0) and (d) Non-Boussinesq (Q = 4.5) 
conditions. In all cases, bright regions in the images indicate the hot upflows and dark regions 
indicate cold downflows in the convective patterns. 



nents due to camera spatial noise. The median value of intensity for all pixels in the image 
is then determined and used as a typical threshold to characterize each pixel as describing 
either hot upflow or cold downflow in the convection pattern. The resulting time series of 
thresholded 1-bit images are used for computing homology. 



III. NUMERICAL SIMULATIONS 



Our direct numerical simulations of the Navier-Stokes equations employ a pseudospec- 
tral code developed by Pesch and co-workers ([3, Q)- The code uses Fourier modes in 
the horizontal direction and appropriate combinations of trigonometric and Chandrasekhar 
functions the satisfy the top and bottom boundary conditions in the vertical direction 10]. 



All runs are performed with six vertical modes and 128 x 128 horizontal Fourier modes in 
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a square domain with side length equal to 16 times the pattern wavelength at convective 
onset. The time step is typically r^/500. For our analysis, the flows are represented by 128 
X 128 images [Figured] (c,d)] of the temperature field or the vertical velocity component. 
The images are typically stored every 2r v . The median value of the flow field quantity 
(temperature or vertical velocity) for each image is determined and used as a threshold to 
characterize each gridpoint as describing hot upflow or cold downflow. Thus, as in the lab 
experiment, the resulting time series of thresholded 1-bit images are used for computing 
homology. 

In the simulations we describe non-Boussinesq effects arising from the temperature de- 
pendence of material properties by a Taylor expansion truncated at leading-order beyond 
the Boussinesq approximation. The simulations are performed at constant mean tempera- 
ture (Th + T r )/2; the expansion is carried out about the mean temperature. In this case, the 
parameter Q ([lOj]) is given by 



where the quantities V% are linear functions of Pr _1 , and the non-Boussinesq coefficients 
7? give the difference of the respective fluid properties across the layer at threshold. In 
non-Boussinesq simulations all the 7? are retained, while in the Boussinesq simulations 7? 
are set to 0. In all simulations, we fix e = 1.4 and set Pr = 0.8. 

IV. RESULTS 

Formally, homology is computed for a topological space X of N dimensions by assigning 
systematically a sequence of abelian groups H k (X) (k = 0, 1,2, ...,N — 1) to X. For our 
purposes, it is sufficient to take H k to be products of the integers: Hk{X) = Z^ x \ where 
the integer dimensions of the groups f3k{X) > are also known as the Betti numbers. 
In this work we focus solely on f3k{X) as the output of the homological analysis; each 
f3k(X) describes a topological property of X. Thus, the net result of this analysis is the 
characterization of X by a set of N non-negative integers. 

From each 1-bit image in time series from either experiments or simulations, two distinct, 
topological spaces are obtained: Xh, where the hot upflow pixels have non-zero value and X c 
where the cold downflow pixels have non-zero value. X c and X^, which are two-dimensional, 
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FIG. 2: Computation of the homology for the experimental data in Figure [U (b) yields a sequence 
of Betti numbers /3, which can be readily interpreted visually. The number of distinct components 
is indicated by the zeroth Betti number for hot upflows (a) (3qh = 65 and for cold downflows (b) 
@0c = 34. (Different colors are used in (a) and (b) to distinguish a given component from its nearest 
neighbors.) The number of holes is given by the first Betti number for hot upflows (c) (5\h = 9 and 
for cold downflows (d) f3\ c = 34. (Each hole is encircled by a red boundary in (c) and (d).) 

are input, in turn, into the homology codes, which subsequently output two Betti numbers 
for each space: fioh, Pih for the hot upflows and j3o c , Pic for the cold downflows. /?oh (A)c) 
counts the number of distinct components, i.e., the number of regions of hot upflow (cold 
downflow) that are separated from similar regions in a given pattern (X c ) [Figure [2] 
(a,b)]. f3\h {Pic) counts the number of holes in the hot upflows (cold downflows) in a given 
pattern Xh (X c ) [Figure [2] (c,d)]. With the package CHomP, computing the homology of X c 
and X h corresponding to each image takes about 1 second. 

Figure [2] shows a striking result: in the experiments, hot upflows are topologically quite 
distinct from cold downflows. Specifically, there are more hot upflow components than 
cold downflow components (/?oft > Aj c )- Moreover, the cold downflow regions contain more 
holes than the hot upflows {(5\ c > /3±h)- This distinction is not revealed using standard 
statistical measures of the pattern. For example, the mean area occupied by upflow is equal 
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FIG. 3: Time series of (a) the zeroth Betti numbers /?ofi, (open circles), (3q c (filled circles) and (b) 
the first Betti numbers f)\h (open diamonds) and I3\ c (filled diamonds) are obtained from laboratory- 
experiments at e = 2.0. Time is scaled by t v ; the time interval between samples is t v /2 



to that occupied by downflow by construction (when the threshold is set to the median pixel 
intensity in the original image.) Wavenumber distributions obtained from Fourier analysis 
of Xh and X c show no discernible differences. 

These measurements of topology are robust to variations in the choice of threshold. The 
choice of the median pixel intensity as the threshold to separate upflows from downflows is 
physically well-motivated but somewhat arbitrary. In practice, any reasonable choice yields 
similar results. For example, for Xh and X c in Figure 2, choosing the mean pixel intensity 
(which is larger than the median intensity by approximately 5% of full scale) as the threshold 
yields nearly identical Betti numbers: j3 h = 65, j3 0c = 34, fi lh = 9, (3\ c = 33. 

Time series of the Betti numbers exhibit fluctuations about well-defined time average val- 
ues (Figure [3]). The fluctuations are primarily a global signature of the complex spatiotem- 
poral behavior of spiral defect chaos. Mean flow induced by curvature in the roll pattern 
leads to regions of local compression or dilatation throughout the pattern. Compression 
often leads to mergering of neighboring rolls while the dilatation results in the formation of 
a new rolls in the pattern; theseprocesses are closely related to secondary instability mech- 



anisms for ideal straight rolls 



10]]. These local events drive further changes in pattern 



curvature, thereby leading to a continually evolving pattern with fluctuating topology. The 
Betti numbers are a global measure of the topological changes, and therefore, are dependent 
on the local processes, for which theories of defect dynamics have been proposed [ 14j]. How 
Betti numbers are related to defect dynamics remains an open question; for our purposes 
here, we focus on the time average values of the Betti numbers (P h,Pih,Poc,Pic), which we 
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FIG. 4: The mean Betti numbers are plotted as a function of the effective optical distance z\ of 
the shadowgraph system in laboratory experiments performed at e = 2. For each data point, the 
median pixel intensity of the raw shadowgraph images was used as the threshold for the homology 
analysis. 



find to be stationary for fixed e. 

The measurements of j3 are robust with respect to flow visualization conditions. It is 
well-known that shadowgraphy can introduce significant nonlinearities and image artifacts 
(e.g., caustics); the strength of these effects depend on the effective optical distance Z\ of 



the shadowgraph system 



11]]. We have checked for possible sensitivity to shadowgraphy 



visualization by conducting a series of experiments where the conditions of the convective 
flow were fixed and image time series were captured for for different values of z\ . Figure H] 
shows that the mean Betti number change only slightly as Z\ is varied over nearly an order 
of magnitude. Additional experimental data (not shown) demonstrate that a change of sign 
in z\ (which changes hot upflows (cold downflows) from bright (dark) to dark (bright)) does 
not affect the determination of /3. 

The differences between the mean Betti numbers for hot upflows j3oh, ft\h an d for cold down 
flows Po c ,Pi c become more substantial as e increased above convective onset (Figure [5]). For 
e < 0.7, the number of upflow components exceeds the number of downflow components 
by a small but statistically significant amount while the number of holes in both types of 
flow patterns are effectively zero. (The onset of spiral defect chaos occurs approximately at 
e = 0.7 in our experiment.) For e > 0.7, the difference in the average component number 
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FIG. 5: The mean zeroth Betti numbers /3o/i (open circles), /?oc (filled circles) and first Betti 
numbers (3\h (open diamonds) and (3\ c (filled diamonds) are shown as a function of e for data from 
laboratory experiments. Each data point is obtained by averaging the Betti numbers from analysis 
of 18000 images corresponding to an observation time of approximately 1800 t v 

grows significantly; at first primarily because the number of hot components grows rapidly 
for 0.7 < e < 2.5 and then because the number of cold components decreases for 2.5 < e < 4. 
The behavior in the number of holes is somewhat different. In the range 0.7 < e < 2.5 the 
number of holes increases significantly in the cold downflows, but only weakly for the hot 
upflows. 

These experimental observations, taken as a whole, suggest the observed asymmetries in 
the Betti numbers may be due to the breakdown of the Boussinesq approximation (and cor- 
responding breaking of the Boussinesq symmetry). To check this hypothesis, we conducted 
two simulations under identical conditions except one simulation was Boussinesq (Q = 0) 
and the other simulation was non-Boussinesq (Q = 4.5). The Betti number time series for 
the Boussinesq simulation shows little distinction in the Betti numbers [Figure [6] (a, b)]. By 
contrast, examination of the same field variable in non-Boussinesq simulation [Figure [6] (c, 
d)] shows distinct differences in the Betti numbers that are qualitatively in agreement with 
the experimental observations. 

Different flow fields extracted from the simulation exhibit similar qualitative behavior in 
the time-averaged Betti numbers [Figure [6] (c, d, e, f)]. We examined both the temperature 
field and the vertical velocity components sampled at the vertical positions z = (the mid- 
plane of the convection cell), z = —0.25 and z = +0.25. (The top and bottom boundaries 
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FIG. 6: Time series of the zeroth Betti numbers Poh (open circles), Pq c (filled circles) and the 
first Betti numbers (3\h (open diamonds) and P\ c (filled diamonds) are obtained from numerical 
simulations at e = 1.4. The midplane temperature field from Boussinesq simulations is used to 
obtain time series of (a) Pq^, Pq c and (b) Pih, Pi c - The midplane temperature field from non- 
Boussinesq simulations is used to obtain time series of (c) Poh, Poc and (d) (3±h, P\ c - The vertical 
velocity component at z = —0.25 from non-Boussinesq simulations is used to obtain time series 
of (e) Poh, Poc and (f) (3ih, Pic- (The midplane is located at z = and the bottom boundary is 
located at z = —0.5.) Time is scaled by t v . 

are located at z — +0.5 and z = —0.5, respectively.) The quantitative values for the 
time- averaged Betti numbers differ weakly between different projections of the convective 
flow. However, every projections exhibit the same qualitative result, namely, in a given 
projection of Boussinesq convection, the time averaged Betti numbers for hot upflows and 
cold downflows are the same whilst each projection of non-Boussinesq exhibits the same 
Betti number asymmetries. 
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V. CONCLUSIONS 



We conclude that the breakdown of the Boussinesq approximation can be readily observed 
in data from convection experiments and simulations by analyzing the topology using com- 
putational homology. It might be argued the use of homology constitutes an "excessive 
use of force" for the two-dimensional patterns analyzed here since the counting of features 
such as components and holes could be accomplished by other means. Nevertheless, the 
simple fact remains that upflow/downflow asymmetries, which had remained unnoticed de- 
spite decades of study of convective flows, were uncovered precisely because of the solid 
mathematical foundation that the homology formalism provides. Moreover, the homology 
analysis outlined here can be readily extended to higher dimensions where less sophisticated 
approaches will likely fail. For example, three-dimensional complexes can be formed from 
the image data used here by creating "time-blocks" of data with two spatial dimensions and 
one time dimension; such data is expected to contain new topological features that capture 
dynamical information ( 151]) . 

Describing the physical mechanism that connects the observed topological asymmetries 
to the breakdown of the Boussinesq approximation remains an open theoretical question. 
There exist physical mechanisms that connect the details of non-Boussinesq effects to the 
hexagonal flow states observed near onset; in particular, the hexagons observed in gases are 
" down- hexagons" , i.e., they have cold downflow in the center of each hexagon. The physical 
mechanism behind this pattern selection is attributed primarily to the fact that as the 
temperature increases in gases, the kinematic viscosity typically decreases and the thermal 
expansion coefficient usually increases [[h]]]- The zeroth Betti numbers of a pattern of M 
ordered hexagons can be easily determined. The centers of each hexagon will be isolated for 
all other cold downflows, yielding f3 0c = M, while the hot upflows around the edges of all 
hexagons will be connected, yielding j3 0h = 1. This behavior of the Betti numbers /5 0c > Poh 
is distinctly different from the results from experiments and simulations at higher e, namely 
floh > A)c- A convincing physical explanation must account for these differences. 

Our results suggest that computational homology might be a useful tool in a wide variety 
of cases in fluid dynamics. For example, in the atmospheric sciences, where extensive use 
is made of the Boussinesq approximation, homological analysis may provide new ways to 
characterize atmospheric data. The use of homology need not be limited to convection; this 
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approach may be applied in any fluid flow where quantitative characterization of complex 
data is needed. 
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FIGURE CAPTIONS 

Figure 1. Images of spiral defect chaos convection are shown from laboratory experi- 
ments (a) & (b) and numerical simulations (c) & (d). Shadowgraph images from the ex- 
periments illustrate the convective flows at (a) e = 1.0 and (b) e = 2.0. The mid-plane 
temperature field is shown at e = 1.4 for simulations carried out under (c) Boussinesq 
(Q = 0) and (d) Non-Boussinesq (Q = 4.5) conditions. In all cases, bright regions in the 
images indicate the hot upflows and dark regions indicate cold downflows in the convective 
patterns. 

Figure 2. Computation of the homology for the experimental data in Figure [I] (b) yields 
a sequence of Betti numbers (3, which can be readily interpreted visually. The number of 
distinct components is indicated by the zeroth Betti number for hot upflows (a) (3^ = 65 
and for cold downflows (b) /3q c = 34. (Different colors are used in (a) and (b) to distinguish 
a given component from its nearest neighbors.) The number of holes is given by the first 
Betti number for hot upflows (c) f3ih = 9 and for cold downflows (d) f3\ c = 34. (Each hole 
is encircled by a red boundary in (c) and (d).) 

Figure 3. Time series of (a) the zeroth Betti numbers floh (open circles), /3q c (filled 
circles) and (b) the first Betti numbers (open diamonds) and /3i c (filled diamonds) are 
obtained from laboratory experiments at e = 2.0. Time is scaled by t v ; the time interval 
between samples is t v /2 

Figure 4. The mean Betti numbers are plotted as a function of the effective optical 
distance z\ of the shadowgraph system in laboratory experiments performed at e = 2. For 
each data point, the median pixel intensity of the raw shadowgraph images was used as the 
threshold for the homology analysis. 
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FIGURE CAPTIONS (cont.) 

Figure 5. The mean zeroth Betti numbers j3 h (open circles), (3 0c (filled circles) and first 
Betti numbers (3^ (open diamonds) and /3i c (filled diamonds) are shown as a function of e for 
data from laboratory experiments. Each data point is obtained by averaging the Betti num- 
bers from analysis of 18000 images corresponding to an observation time of approximately 
1800 t v . 

Figure 6. Time series of the zeroth Betti numbers (3 h (open circles), (3 0c (filled circles) 
and the first Betti numbers fi\h (open diamonds) and fi\ c (filled diamonds) are obtained 
from numerical simulations at e = 1.4. The midplane temperature field from Boussinesq 
simulations is used to obtain time series of (a) (3^, (3q c and (b) /3ih, f3\ c . The midplane 
temperature field from non-Boussinesq simulations is used to obtain time series of (c) floh, 
Poc and (d) /3ih, (3i c - The vertical velocity component at z = —0.25 from non-Boussinesq 
simulations is used to obtain time series of (e) /3 /i, Po c and (f) (3^, fi lc . (The midplane is 
located at z = and the bottom boundary is located at z — —0.5.) Time is scaled by t v . 
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Figure 2. 
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Figure 3. 
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Figure 5. 
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Figure 6. 
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